Bifurcation and hybrid control of a discrete eco-epidemiological model with Holling type-III

In this paper, a three dimensional discrete eco-epidemiological model with Holling type-III functional response is proposed. Boundedness of the solutions of the system is analyzed. Existence condition and stability of all fixed points are discussed for the proposed model. Furthermore, we obtained the transcritical bifurcation surfaces of the system by bifurcation theory. Based on the explicit criteria for the Neimark Sacker bifurcation and flip bifurcation, we obtained that the system undergoes these two types of bifurcations at the positive fixed point. Then we apply a hybrid control strategy that based on both parameter perturbation and a state feedback strategy to control the Neimark-Sacker bifurcation. Finally, some numerical simulations are carried out to support the analytical results.


Introduction
Eco-epidemiology is considered to be a relatively new branch of mathematical biology that studies both ecological and epidemiological issues simultaneously.In 1986, Anderson and May [1] firstly considered disease factor in the predator-prey system.They coupled the epidemiological model developed by Kermack and Mckendrick [2] to a Lotka-Volterra predatorprey model [3,4].Chattopadhyay and Arino [5] coined the term 'eco-epidemiology' [6,7], and proposed a predator-prey epidemiological model with disease spreading in prey, and one of its simplified models takes the form as follows where S, I, Y represent the populations of susceptible prey species, infected prey species and predator species, respectively.In 2001, Chattopadhyay et al [8] studied the practical problem of Pelicans at risk in the Salton Sea based on this model.The Salton Sea, the largest inland water body in California, USA, is a eutrophic salty lake.In the summer, the weather reaches 128 degrees Fahrenheit and the water evaporates very quickly [9].This process increases the salinity of the Salton Sea and reduces oxygen levels, as saltwater is more difficult to combine with oxygen than freshwater.Four types of fishes are very common in Salton Sea, namely Tilapia, Corvina, Croaker, and Sargo.Among them, Tilapia is the most abundant because of its amazing reproductive rate [10].The Salton Sea is the main habitat for many migratory birds such as pelicans, but thousands of water birds (most of which are pelicans) and fish have died.The real cause of this is not yet clear, but there is growing evidence pointing to toxic algal blooms.Algal blooms grow and die very quickly, and in doing so, it draws oxygen from seawater.A lack of oxygen in the tissues of infected fish can lead to outbreaks of botulism, and the later stages of the disease cause the infected fish to rise closer to the surface in search of oxygen.As reported by US Geological Survey National Wildlife Health Center, many white pelicans died from botulism type C between 1978 and 2003 [11].In 1996, over 8500 white Pelicans died due to infection with type C botulism in the Salton Sea.Millions of dead or sick fish can transmit botulism to Pelicans and other migratory birds that feed on fish.As pelicans prey on vulnerable fish, the ingestion of botulinum has lead to the development of avian botulism [12].In the past, many research has been done considering various aspects of interaction and interrelation of Tilapia fish, botulism type C, and Pelican birds.System (1) and its extensions have been investigated under different conditions and functional response.Such as the eco-epidemiological predator-prey model with diseases in the prey [13][14][15][16][17], with a disease in the predator [18][19][20][21], with disease in both populations [22][23][24], with delay [25,26] and with Holling type functional response [27][28][29][30].
However, these studies mostly focus on continuous time system, and rarely involve discrete time system.It is also important to consider discrete-time models.Firstly, discretization of continuous time model is the basis of obtaining numerical approximate solution [31].Secondly, due to the fact that statistical data of epidemics are collected in discrete time, it is more convenient and accurate to describe epidemics with a discrete-time model [32,33].Thirdly, many species, such as monocarpic plants, and semelparous animals have independent and non-overlapping generations, and their births occur during the regular breeding season.Their interactions are described using difference equations or discretetime models [34].Moreover, even a single-species discrete-time model can exhibit bifurcation, chaos, and more complex dynamic behaviors.Recently, discrete-time models have received more and more attention, see [35][36][37][38][39][40] and the references cited therein.Whereas, these studies mainly focus on two-dimensional discrete-time systems, with relatively few studies on the dynamics of three-dimensional discrete-time systems, and even fewer studies on the dynamics of three-dimensional discrete-time eco-epidemiological systems.
Lately, several works related to discrete-time eco-epidemiological models have appeared in the literature.In [41], the authors considered the discretization system of system (1) with ratio-dependent Michaelis-Menten functional response.Then, in [42], the authors considered a discrete system with saturated incidence rate based on paper [41], and obtained the stability, bifurcation and chaos of the system.In [43], the authors considered the discretization system of system (1) with Holling type-II functional response.But there has been limited study on three-dimensional discrete-time systems with Holling type-III functional response, due to its highly non-linear nature.This motivated us to consider a discrete-time eco-epidemiological model with Holling type-III functional response incorporating disease in prey to study the interaction between Tilapia and the Pelican.In this article, the existence and stability of fixed points, bifurcation analysis and a hybrid control strategy are discussed.Through this article, we make the following assumptions.
(A1) The disease is spread among the prey population Tilapia only and the disease is not genetically inherited.The infected population does not recover or become immune.The prey is divided into two classes, susceptible x and infective y.The disease incidence rate adopts bilinear incidence rate βxy, and β is the transmission coefficient.
(A2) In the absence of disease, the prey population Tilapia grows according to the logistic law with intrinsic growth rate r and carrying capacity K.Only susceptible prey is capable of reproducing, the infective preys cannot produce offsprings due to the disease or by predation before having the possibility of reproducing.However, as infected prey still consumes resources, it still contributes to carrying capacity.Thus, the growth rate for the prey is rx(1 − (x + y)/K).
(A3) The predator population Pelican consumes only the infected Tilapia.This is conform to the fact that the infected individuals are less active and can be more easily captured [8].
(A4) The predator consumes prey following Holling type-III functional response.In fact, if the predator actively seeks out large concentration of prey the Holling type-III function f(X) = aX 2 /(m + X 2 ) is more appropriate, where a is the maximum predation rate and m is half saturation constant.Since the slope of this function goes to zero for small values of X(f'(0)=0), it may be suspected that the food chain will be destabilized if prey concentration becomes too small [44].However, after a certain value of X, the predators increase their feeding rates until some saturation level is reached (f(X) ! a when X ! 1).Holling type III functional response is mostly used when the number of predator encounters the prey population with a very lower amount due to unavailability of prey but when the prey becomes available then the response behaves like the Holling II type functional response (f(X) = aX/(m + X)) [45].
Although both Holling II and III functional responses are approaching an asymptote, the former is decelerating and the latter is sigmoid [46].Holling [46] suggested that the type II responses are characteristic of predators, which have no learning ability or when give only one type of prey for which to search, whereas type III responses are characteristic of vertebrate predators where switching and learning are more common.We show the Holling type-II and III function in Fig 1.Under the above assumptions, we formulate the following discrete-time eco-epidemiological model where x n , y n , z n � 0 represent the densities of susceptible prey population Tilapia, infected prey population Tilapia and their predator population Pelican at time n, respectively.Here parameters r, K, β, m, a, b, c, d are all nonnegative constants and their biological meanings are given Table 1 The rest of this paper is organized as follows.In Section 2, The boundedness of solutions of system (2) is given.In Section 3, we study the existence and stability of the fixed points.In Section 4, the transcritical bifurcation of boundary fixed points and the flip and Neimark-Sacker bifurcation of positive fixed points are analyzed.In Section 5, we use a hybrid control strategy to control the Neimark-Sacker bifurcation and the flip bifurcation.Finally, numerical simulations and conclusions are given in Section 6.

Boundedness
For system (2) we always assume that all initial values are non-negative and all the parameters are positive.The boundedness of solutions of system (2) is given by the following lemma.
Lemma 1 Let μ = min(c, d), then all the solutions of system (2) will lie in the region for all initial values x 0 , y 0 , z 0 � 0 as n !+ 1. Proof 1 For system (2), we always assume that x 0 , y 0 , z 0 � 0. Because the environmental carrying capacity of the prey population is K, therefore x n � K. Let x n + y n + z n = M n , adding all the equations of system (2), we get where μ = min{c, d}.
If there exists l 2 N such that M l+1 > M l , then we get Hence, we have that M l � M � , where M � = (r + 1)K/μ.We claim that M n � M � for all n � l.Otherwise, we assume that there exists a q 2 N such that M q > M � , where q � l.Let q be the smallest integer such that M q > M * , then M qÀ 1 � M * .From inequality (3) we get M q � M * , which is a contradiction.Hence the claim is hold. If Otherwise, we assume that � M > M * .Taking limit of equation for large enough n 2 N, which is a contradiction.Hence the claim is proved.Thus, all the solutions of system (2) will lie in the region B.

Existence of fixed points and stability
In this section, we give the sufficient and necessary conditions for the existence of fixed points for system (2), and analyze their properties.Defined R 0 = βK/c.Because β is the transmission coefficient, so βK represents the number of newly infected individuals when all the prey populations are susceptible at the beginning of the disease.For c is the death rate of infected prey because of natural death and disease induced mortality, so 1/c is the duration of infection of an infected prey.Therefore, R 0 is the disease reproduction number in the prey.
Firstly, on the sufficient and necessary conditions for the existence of the nonnegative fixed points of system (2), we give the the following theorem.
Since the second equation in (7) implies y n = 0 if x n = 0, the third equation in ( 7) implies z n = 0 if y n = 0. Therefore, for the boundary equilibrium point, we only need to consider the case of x n 6 ¼ 0, y n 6 ¼ 0, z n = 0.If x n 6 ¼ 0, y n 6 ¼ 0, z n = 0 and R 0 � 1, the second equation in (7) holds only if y n = 0. Thus, there are only two solutions (0, 0, 0) and (K, 0, 0) of ( 7) when R 0 � 1.If x n 6 ¼ 0, y n 6 ¼ 0, z n = 0 and R 0 > 1, (7) has a solution ðx; ŷ; 0Þ, where x; ŷ are given in (5). When from the third equation of ( 7) we can see that when ab � d, the third equation of ( 7) holds only if z n = 0.When ab > d, if z n 6 ¼ 0, from the third equation of ( 7), we can get y n ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dm=ðab À dÞ p . From the first equation of ( 7), we can get x n ¼ K À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dm=ðab À dÞ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dm=ðab À dÞ the second equation of ( 7) not holds.
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dm=ðab À dÞ p bðr þ bKÞ=ðcrÞ, the second equation of ( 7) holds only if z n = 0.If z n = 0, then x n = 0, y n = 0 or x n = K, y n = 0 or x n ¼ x, y n ¼ ŷ.Thus, there are only three solutions (0, 0, 0), (K, 0, 0) and ðx; ŷ; 0Þ of ( 7 When ab > d and R 0 > 1 þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dm=ðab À dÞ p bðr þ bKÞ=ðcrÞ, if z n 6 ¼ 0, from the third equation of (7), we can get y n ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dm=ðab À dÞ p > 0. From the first equation of ( 7), we get From the second equation of ( 7), we get z n ¼ bmðbx n À cÞ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dmðab À dÞ p ¼ bcm ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dmðab À dÞ so ( 7) has a unique positive solution (x * , y * , z * ), where x * , y * , z * are given in (6).If z n = 0, then x n = 0, y n = 0 or x n = K, y n = 0 or x n ¼ x, y n ¼ ŷ.Thus, there are exactly four fixed points E 0 (0, 0, 0), E 1 (K, 0, 0), E 2 ðx; ŷ; 0Þ and For a discrete dynamical system on R 3 , let the Jacobian matrix of this system at a fixed point (x, y, z) be J(x, y, z).We denote the three eigenvalues of J(x, y, z) by λ i (i = 1, 2, 3).We recall some concepts for a discrete dynamical system on R 3 [47].If |λ i | < 1 for all eigenvalues, then (x, y, z) is called a sink and is locally asymptotically stable; if |λ i | > 1 for some eigenvalues and |λ j | < 1 for the others, then (x, y, z) is called a saddle and is unstable; if |λ i | > 1 for all eigenvalues, then (x, y, z) is called a source and is unstable; if |λ i | = 1 for any eigenvalue, then (x, y, z) is called non-hyperbolic.
Theorem 2 (1) The fixed point E 0 is always unstable.
(2) For the fixed point E 1 , we have the following conclusions.
Proof 3 Because the Jacobian matrix of system (2) at E 0 is given by one of the eigenvalue of matrix J(E 0 ) is λ 1 = 1 + r > 1, so the fixed point E 0 is always unstable.The Jacobian matrix of system (2) at E 1 is given by Thus, E 1 is non-hyperbolic in this case.
Lemma 3 (Jury-criterion [49]) For the equation λ 3 + a 2 λ 2 + a 1 λ + a 0 = 0, all roots lie within the unit disk if and only if the following conditions are satisfied, where a 0 , a 1 , a 2 are real numbers.
In the follows, we use Lemma 3 to analyze the stability of the positive fixed point E 3 .Theorem 4 When E 3 exists, E 3 is locally asymptotically stable if and only if the following inequations hold true Proof 5 The Jacobian matrix of system (2) at E 3 is The corresponding characteristic equation of J(E 3 ) is where b i (i = 0, 1, 2) are defined in (9).According to Lemma 3, we get that E 3 is locally asymptotically stable if and only if condition ( 8) is satisfied.

Bifurcations analysis
In this section, we investigate possible bifurcation of system (2).

Transcritical bifurcation
In this subsection, we will give the transcritical bifurcation of system (2).Proof 6 From Theorem 2, we know that the eigenvalues for E 1 are 1 − r, βK − c + 1 and 1 − d.In O 1 , we require r < 2 to ensure |1 − r| < 1.In this case, the unique non-hyperbolic case is exact βK − c + 1 = 1, which corresponds to transcritical bifurcation surface O 1 .For the critical case βK − c = 0, that is R 0 = 1, system (2) has only two fixed points E 0 and E 1 .It is not hard to check that bK À c À bðr þ bKÞ r ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi which means that the parameter goes into Λ 1 [ Λ 2 when R 0 increases from 1. Thus, the third fixed point E 2 appears by Theorem 1, i.e., transcritical bifurcation happens for E 1 .We observe that E 2 is sufficiently close to E 1 when 0 < R 0 − 1 � 1.
When E 2 exists, from Theorem 3, we know that one of the eigenvalues for E 1 is w 1 .If F E 2 ðÀ 1Þ > 0 and C E 2 < 1, i.e., 2rc/(βK) − 4 < r(βK − c) < rc/(βK), r < 4βK/c, the modulus of the other two eigenvalues with E 2 is less than 1.In this case, the unique non-hyperbolic case is exact w 1 = 1, which corresponds to transcritical bifurcation surface O 2 .When the parameter crosses O 2 into Λ 3 , the forth fixed point E 3 appears by Theorem 1, i.e., transcritical bifurcation happens for E 2 .We observe that E 3 is sufficiently close to E 2 when 0 < R 0 À ð1 þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dm=ðab À dÞ is the characteristic polynomial of J(E 3 ).From Theorem 1, we know that when E 3 exists, we have ab − d > 0 and βx * − c > 0. Thus, we can judge that This indicates that a fold bifurcation not happens for E 3 of system (2).Thus, in the following, we will investigate the flip bifurcation and Neimark-Sacker bifurcation for E 3 of system (2).

Flip bifurcation
In this subsection, we will discuss parametric conditions under which the unique positive fixed point of system (2) undergoes a flip bifurcations.For this purpose, an explicit criterion for flip bifurcation is implemented without finding the eigenvalues of the system.The criterion is formulated using a set of simple equalities or inequalities that consist of the coefficients of the characteristic equation derived from the Jacobian matrix.Next, let's introduce this criterion first [50].Consider an n-dimensional map x n+1 = f μ (x n ), where x n+1 , x n 2 R n and μ 2 R is a parameter.Assume that f has a fixed point x 0 and the characteristic polynomial of an n-dimensional map f μ at x 0 takes the form where σ 0 = 1 and σ i = σ i (μ, v) (i = 1, 2, � � �, n), μ is the bifurcation parameter, and v is the control parameter or the other to be determined.Consider the sequence of determinants D � 0 ðm; vÞ ¼ 1, D � 1 ðm; vÞ; � � � ; D � n ðm; vÞ, where Lemma 4 [50] Assume that f μ has a fixed point x 0 .A flip bifurcation takes place at μ = μ 0 if and only if the following conditions (H1) Eigenvalue assignment: P m 0 ðÀ 1Þ ¼ 0, D � nÀ 1 ðm 0 ; vÞ > 0, P m 0 ð1Þ > 0 and D � j ðm 0 ; vÞ > 0 for j = n − 2, n − 4, � � �, 1 (resp.2), when n is odd (resp.even).(H2) Transversality condition: are satisfied, where s 0 i stands for the derivative of σ i (μ) with respect to μ at μ = μ 0 .When n = 3, and we choose r as the perturbation parameter, in the following theorem, we give the parametric conditions for the flip bifurcation takes place at r = r 0 for E 3 of system (2).
Theorem 6 The fixed point E 3 (x * , y * , z * ) of system (2) undergoes a flip bifurcation at r = r 0 if the conditions are satisfied, where b i (i = 0, 1, 2) are defined in ( 9), b 0 i is derivative of b i (r) with respect to r at r = r 0 , and r 0 is a possible real root of equation 1 − b 2 (r) + b 1 (r) − b 0 (r) = 0.
Proof 7 For n = 3 and r 0 is the perturbation parameter.According to the criterion introduced in Lemma 4, if the conditions (H1) and (H2) are satisfied, then a flip bifurcation occurs at r 0 for system (2).That is Then we get the conditions (10).

Neimark-Sacker bifurcation
In this subsection, we discuss parametric conditions under which the unique positive fixed point E 3 of system (2) undergoes a Neimark-Sacker bifurcations.For this purpose, an explicit criterion for Neimark-Sacker bifurcation is implemented without finding the eigenvalues of the system.We state the explicit criterion as follow.
When n = 3 and r is taken as the bifurcation parameter, we give the parametric conditions for the Neimark-Sacker bifurcation takes place at r = r 0 for E 3 of system (2) in the following.
Theorem 7 The fixed point E 3 (x * , y * , z * ) of system (2) undergoes a Neimark-Sacker bifurcation at r = r 0 if the conditions are satisfied, where b i (i = 0, 1, 2) are defined in (9), and r 0 is a possible real root of equation Proof 8 For n = 3 and r is the perturbation parameter.According to the criterion introduced in Lemma 5, if the conditions (C1), (C2) and (C3) are satisfied, then Neimark-Sacker bifurcation occurs at r 0 for system (2).That is Then we get the conditions (11).

Bifurcation control
In order to prevent the serious damage or even extinction of the population caused by infectious diseases, a stable positive fixed point may be needed to maintain the sustainable development of the eco-epidemiological system.That is, it is better for the positive fixed point E 3 to be an asymptotically stable of system (2) if it exists.Therefore, we would like to take certain control measures to avoid the happening of bifurcations.For this purpose, in this section we provide for system (2) a hybrid control, which is based on feedback control strategy and parameter perturbation (see [52,53]).Corresponding to the system (2), we construct a controlled system where 0 < θ < 1, and θ is called a control parameter.Obviously, controlled system (12) is exactly (2) if θ = 1.System (12) has same fixed points as (2).The Jacobian matrix for ( 12) at The corresponding characteristic equation of J con is where Then by Lemma 3 we have the following theorem.Theorem 8 When E 3 exists, the unique positive fixed point E 3 of system ( 12) is locally asymptotically stable if and only if where c i (i = 0, 1, 2) are defined in (14).Proof 9 Clearly, we only need to prove that the modulus of all eigenvalue is less than 1 for characteristic equation (13), which is equivalent to by the conditions given in Lemma 3, then we get the conclusion in this theorem.

Numerical simulations and conclusions
In this section, we give the phase portraits, bifurcation diagrams, and Lyapunov exponents to illustrate our theoretical results numerically.We selected parameters based on the reference [8,43] that studied the interaction between Pelican and Tilapia in the Salton Sea, and their biological meanings are stated in Table 1.
In  corresponding to the Neimark-Sacker bifurcation point.After that, we observe that some Lyapunov exponents are positive and some are negative, confirming the existence of the chaotic regions in the parametric space.
In  (2).In this case, from the condition (15), we get Inequ1 : j À 3:95397y 2 þ 7:95999y À 4:04529j < 3:95397y 2 À 7:88133y þ 4:04529; Inequ2 : j15:58y 3 À 23:75y 2 þ 8:30y À 0:09j < j À 15:63y 4 þ 31:79y 3 À 24:42y 2 þ 8:40y À 0:09j; Inequ3 : j À 3:95397y 2 þ 4:01933y À 1:04529j < 1: According to Theorem 8, we know that when condition (15) is satisfied, we can control E 3 to be a stable fixed point.After a simple calculation, we obtained that E 3 is stable when 0.0205565 < θ < 1.Thus, we take θ = 0.99999, and give the time-series graph in Fig 8(iv)-8(vi) for system (12).We observe that the Neimark-Sacker bifurcation of system (2) is controlled effectively.We give the stable region of system (12) for β, θ 2 (0, 1) in Fig 9.In this paper, we have constructed the mathematical model to describe an interaction between Tilapia as a prey, Pelicans as predator, and botulism type C as the cause of disease in Tilapia.We discuss the dynamical behaviors of the system (2), establish that the system solution is bounded and get that the system has at most four fixed point depending on values of the system parameters.We obtained a threshold value R 0 = βK/c such that R 0 � 1 leads to the complete disappearance of infected prey from the ecosystem and the eradication of the disease from the prey population.If the death rate of infected prey c is high and transmission rate from uninfected prey to infected prey β is low, then the chance of disease eradication from the ecosystem will be high.This means that we can control the spread of the disease by killing botulism type C with drugs.The local asymptotic stability of different fixed points are discussed here.The fixed point E 0 (0, 0, 0) is always unstable, which implies that the system can never be collapsed for any values of the system parameters.The fixed point E 1 (K, 0, 0) and E 2 ðx; ŷ; 0Þ are locally asymptotically stable under some parametric restrictions.There exists a set of values of the system parameters for which the positive fixed point E 3 (x * , y * , z * ) is locally asymptotically stable, i.e. both the populations can survive with positive density level.We notice that the stability of the fixed points is greatly influenced by the model parameters.We show that under certain parametric conditions system (2) undergoes transcritical bifurcations at boundary fixed points E 1 and E 2 , using bifurcation theory.Further, by the explicit criteria for a flip bifurcation and a Neimark-Sacker bifurcation, we prove that the system undergoes both flip and Neimark-Sacker bifurcations at the fixed point E 3 under some parametric conditions.From the ecological point of view, flip bifurcation is associated with the emergence of chaotic behavior, demonstrating the evolution of the prey and predator populations.An invariant curve bifurcates from the fixed point, meaning that predator and prey can coexist in a stable way and reproduce their densities.The dynamics on the invariant curve may be either periodic or quasi-periodic [40].The Maximum Lyapunov exponents are numerically computed to confirm further the complexity of the dynamical behavior.Finally, we use the hybrid method to control the Neimark-Sacker bifurcation at fixed point E 3 .The result indicate that the nonlinear dynamics of such eco-epidemiology model not only depend on more bifurcation parameters but also are very sensitive to parameter perturbations, which are important for the control of biological species or infectious diseases.Finally, numerical simulations are carried out to confirm the validity of the theories and the effectiveness of the control method.

Fig 7 ,
we give the Neimark-Sacker bifurcation diagrams for system (2) in the (r, x n )plane, (r, y n )-plane and (r, z n )-plane.The Maximum Lyapunov exponents corresponding to Fig 7(i)-7(iii) are calculated and plotted in Fig 7(iv), which are consistent with the bifurcation diagram.When r 2 (1.7, 1.791065), the Maximum Lyapunov exponents are negative, which denote the system is stable.At r = 1.791065, the Maximum Lyapunov exponent is 0, which is

Table 1 . Description and estimation of parameters in system (2).
for E 1 and E 2 , respectively.That is, transcritical bifurcation happens for E 1 when the parameter crosses O 1 into Λ 1 [ Λ 2 and E 2 appears near E 1 ; transcritical bifurcation happens for E 2 when the parameter crosses O 2 into Λ 3 and E 3 appears near E 2 .